Jump to Empirical Spatial Indicator Analysis
This vignette provides an overview of spatial indicators used to estimate the location, range, occupancy, and aggregation of populations in the context of fisheries management.
Spatial indicators are configured to work with DATRAS survey data. Load in some survey data and process it first before parsing through spatial indicator functions.
Load DATRAS exchange data
Here we load NS-IBTS Quarter 1 survey data for years 2018 to 2022.
yrs <- c(2000:2022)
qrs <- c(1)
srv <- "NS-IBTS"
hh <- icesDatras::getDATRAS(record = "HH", srv, years = yrs, quarters = qrs)
hl <- icesDatras::getDATRAS(record = "HL", srv, years = yrs, quarters = qrs)
ca <- icesDatras::getDATRAS(record = "CA", srv, years = yrs, quarters = qrs) # for completeness, but not used
Remove duplicate rows in exchange data.
hh <- unique(hh)
hl <- unique(hl)
ca <- unique(ca)
ICES Statistical Rectangles
In the exchange data (in hh and hl) we know
the ICES rectangle that each haul was conducted in, but we do not have
the associated division. To be able to subset data to a specific stock
region, we need to download ICES Statistical
Rectangle shapefile (Quick Downloads > ICES StatRec mapped to ICES
Area) and append information on ICES divisions to our exchange data.
A download is available on GitHub and has been converted to a
.rds file which is loaded in here.
load(url("https://github.com/peterjosephkidd/SpatIndAssess/raw/main/Data/ICES%20Rect/ices_rect.rds"))
Add ICES Divisions
area_div <- dplyr::distinct(ices_rect[c("ICESNAME", "Area_27", "Shape_Area")])
hh <- merge.data.frame(hh, area_div, by.x = "StatRec", by.y = "ICESNAME")
ca <- merge.data.frame(ca, area_div, by.x = "AreaCode", by.y = "ICESNAME")
Basic Data Processing
Remove invalid hauls:
hh <- filter(hh, !HaulVal %in% c("I", "P")) # p = partly valid, it is deprecated
-9 is a placeholder for NAs. Change to NA
for TotalNo column:
hl$TotalNo[hl$TotalNo == -9] <- NA
Create haul.id:
hh$haul.id <- as.character(
paste(hh$Year, hh$Quarter, hh$Country, hh$Ship, hh$Gear, hh$StNo, hh$HaulNo,
sep = ":"))
hl$haul.id <- as.character(
paste(hl$Year, hl$Quarter, hl$Country, hl$Ship, hl$Gear, hl$StNo, hl$HaulNo,
sep = ":"))
ca$haul.id <- as.character(
paste(ca$Year, ca$Quarter, ca$Country, ca$Ship, ca$Gear, ca$StNo, ca$HaulNo,
sep = ":"))
Merge haul information from hh with landings information
in hl.
## Refine columns
m <- hh[c("haul.id", "Year", "Quarter", "Month", "Survey","Country", "Ship",
"Gear", "GearEx", "DoorType", "HaulDur", "HaulNo", "StNo", "SweepLngt",
"StatRec", "Area_27", "ShootLong", "ShootLat", "HaulVal", "Depth")]
## Merge HL with HH
hlhh <- merge(hl, dplyr::distinct(m),
c("haul.id", "Year", "Quarter", "HaulNo", "StNo", "Gear", "GearEx",
"DoorType","Ship", "SweepLngt", "Country", "Survey"))
rm(m)
Now that DATRAS survey data has been loaded and formatted, choose a stock to investigate and filter data accordingly. This example looks at European plaice (Plueronectes platessa) in Subarea 4 (North Sea) and Subdivision 20 (Skagerrak; ple.27.420).
Set parameters
# ple.27.420
species <- c("plaice", "Pleuronectes platessa")
species_aphia <- findAphia(species[2], latin = T)
stk_divs <- c("4.a", "4.b", "4.c", "3.a.20")
Filter data
Basic checks
# Source Spatial Indicator functions from GitHub
i <- 1
source("https://raw.githubusercontent.com/peterjosephkidd/SpatIndAssess/main/Functions/spatinds_funs.R", print.eval = F)
The mean longitudinal and latitudinal location of the population. Changes in CG indicate whether the mean location of the population is shifting eastward/westward (CGx) or northward/southward (CGy).
cg <- coginis(hlhh, yrs, qrs, species_aphia, stk_divs,
iso = F, inertia = F, # toggle off outputting two another spatial indicators for now
density = T) # weight hauls by density of catch at each sample site (F therefore uses binary presence-absence data)
head(cg)
Year Quarter CoG (x) CoG (y)
1 2000 1 6.35 55.8
2 2001 1 4.74 54.7
3 2002 1 4.25 54.7
4 2003 1 4.16 54.6
5 2004 1 4.95 55.1
6 2005 1 5.87 55.8
Visualise changes in cg over time
cg_p <-cogplot(cg,
# set grid = ices_rect for ICES rectangles -- takes time to load
grid = ices_rect, areas = stk_divs,
xlim = c(-2, 6), ylim = c(52, 60))
suppressWarnings(ggplotly(cg_p, tooltip = "text"))
Inertia describes the dispersion/variance of the population around its centre of gravity. High values of inertia indicate that the population is widely spread across space.
inert <- coginis(hlhh, yrs, qrs, species_aphia, stk_divs,
iso = F, inertia = T, # inertia toggled on
density = T) #
head(inert)
Year Quarter CoG (x) CoG (y) Inertia
1 2000 1 6.35 55.8 14.8
2 2001 1 4.74 54.7 12.3
3 2002 1 4.25 54.7 12.8
4 2003 1 4.16 54.6 13.5
5 2004 1 4.95 55.1 19.7
6 2005 1 5.87 55.8 19.6
Plot the trend over time
in_p <- ggplot(data = inert, aes(x = Year, y = Inertia)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Inertia (I) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA))
in_p
A convex hull polygon is drawn around occurrence points (i.e. sample sites with catch of the target species > 0). EOO is the area of the polygon. High EOO indicates that the population is spread over a large geographical area.
eoo <- chullarea(hlhh, yrs, qrs, species_aphia, stk_divs)[1:3] # this function needs tidying up, ignore other columns for now
head(eoo)
Year Quarter convex_hull_area
1 2000 1 89.3
2 2001 1 85.9
3 2002 1 90.0
4 2003 1 87.0
5 2004 1 87.5
6 2005 1 88.3
Plot the trend over time
eoo_p <- ggplot(data = eoo, aes(x = Year, y = convex_hull_area)) +
geom_line() +
geom_point() +
theme_minimal() +
ylab("Extent of Occurence") +
labs(title = "Extent of Occurence (EOO)",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA))
eoo_p
Similar to EOO. But instead of a convex hull, ELA calculates the area of an ellipse that encompasses 95% of the occurrence sites.
print(stk_divs)
[1] "4.a" "4.b" "4.c" "3.a.20"
ela <- ellarea(hlhh, yrs, qrs, species_aphia, stk_divs = stk_divs)
head(ela)
Year Quarter Ellipse Area
1 2000 1 86.2
2 2001 1 75.0
3 2002 1 86.6
4 2003 1 92.0
5 2004 1 102.0
6 2005 1 109.4
ela_p <- ggplot(data = ela, aes(x = Year, y = `Ellipse Area`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Ellipse Area (ELA) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA))
ela_p
Indicates the proportion of area occupied by the target species using binary presence-absence data.
The proportion of the ICES rectangles surveyed with a catch of the target species > 0.
popr <- pa_rect(hlhh, yrs, qrs, species_aphia, stk_divs)
head(popr)
# A tibble: 6 × 5
Year nrects Quarter nrects_p PosAreaR
<dbl> <int> <chr> <dbl> <dbl>
1 2000 163 1 127 0.779
2 2001 165 1 128 0.776
3 2002 162 1 131 0.809
4 2003 165 1 123 0.745
5 2004 165 1 131 0.794
6 2005 166 1 131 0.789
popr_p <- ggplot(data = popr, aes(x = Year, y = `PosAreaR`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Proportion of Presence in ICES Rectangles (POPR) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA)) +
ylab("Proportion of Presence (Rectangle)")
popr_p
The proportion of hauls with catch of the target species > 0.
poph <- pa_haul(hlhh, yrs, qrs, species_aphia, stk_divs)
head(poph)
# A tibble: 6 × 5
Year no_haul.ids Quarter pr_hauls PosAreaH
<dbl> <int> <chr> <dbl> <dbl>
1 2000 367 1 256 0.698
2 2001 411 1 288 0.701
3 2002 401 1 297 0.741
4 2003 398 1 276 0.693
5 2004 355 1 246 0.693
6 2005 371 1 261 0.704
poph_p <- ggplot(data = poph, aes(x = Year, y = `PosAreaH`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Proportion of Presence in Survey Hauls (POPR) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA)) +
ylab("Proportion of Presence (Haul)")
poph_p
The mapdis() function visualises the above range and
occupancy indicators on a single map for one timepoint (e.g. year and
quarter). The function can equally be used to visualise a single
indicator.
mapdis(hlhh, yrs = 2022, qrs, species_aphia, stk_divs, ices_rect, # data & parameters
cog = T, inertia = T, EOO = T, ELA = T, # indicator toggles
density = T, # weight samples
title = "Plaice (Pleuronected platessa)\nNS-IBTS",
xlim = c(-5,11), ylim = c(50, 62)) # plotting window
Derived from a Lorenz curve. Ranges from 0 to 1, with 1 indicating that the population is uniformly distributed across surveyed rectangles, and 0 indicating that the population was recorded in one rectangle.
lordat <- lorenz_data(hlhh, yrs, qrs, species_aphia, stk_divs)
lorenz_plot(lordat) +
theme_minimal() +
theme(panel.border = element_rect(colour = "black", fill = NA)) # function needs tidying
gni <- Gini(lordat)
head(gni)
# A tibble: 6 × 3
Year Quarter `Gini Index`
<dbl> <chr> <dbl>
1 2000 1 0.219
2 2001 1 0.199
3 2002 1 0.231
4 2003 1 0.181
5 2004 1 0.266
6 2005 1 0.256
gni_p <- ggplot(data = gni, aes(x = Year, y = `Gini Index`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Gini Index Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA))
gni_p
Represents the proportion of the population present in 95% of the area. Ranges from 0 to 0.95, with 0 indicating that all individuals were recorded in 5% of the surveyed area (high aggregation) and 0.95 indicating that 95% of the population were recorded in 95% of the rectangles surveyed (uniform distribution).
D <- d95(lordat)
D95 is an estimate of the proportion of the population that exists in 95% of surveyed rectangles.
d95_p <- ggplot(data = D, aes(x = Year, y = `D95`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Proportion of Catch in 95% of ICES Rectangles in Stock Region (D95) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA))
d95_p
Derived from cumulative frequency data. Measures the area occupied by the population, taking into account variations in fish density values. High values of SA indicate homogeneous spatial distribution.
sa_data <- spreadingarea_data(hlhh, yrs, qrs, species_aphia, stk_divs)
sa <- sa_data %>%
group_by(Year) %>%
summarise("Spreading Area" = spreadingarea_calc(TotalNo_Dur)) %>%
mutate(Quarter = paste(as.character(sort(unique(sa_data$Quarter))), collapse = ", ")) %>% # add quarters
relocate(Year, Quarter)
head(sa)
# A tibble: 6 × 3
Year Quarter `Spreading Area`
<int> <chr> <dbl>
1 2000 1 80.4
2 2001 1 87.0
3 2002 1 97.7
4 2003 1 74.4
5 2004 1 95.9
6 2005 1 95.3
sa_p <- ggplot(data = sa, aes(x = Year, y = `Spreading Area`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Spreading Area (SA) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA))
sa_p
The area that would be covered by a population with homogeneously distributed density. It is equal to the mean density per individual.
ea <- sa_data %>%
group_by(Year) %>%
summarise("Equivalent Area" = equivalentarea(TotalNo_Dur)) %>%
mutate(Quarter = paste(as.character(sort(unique(sa_data$Quarter))), collapse = ", ")) %>% # add quarters
relocate(Year, Quarter)
head(ea)
# A tibble: 6 × 3
Year Quarter `Equivalent Area`
<int> <chr> <dbl>
1 2000 1 45.9
2 2001 1 58.5
3 2002 1 73.8
4 2003 1 41.5
5 2004 1 72.5
6 2005 1 78.9
ea_p <- ggplot(data = ea, aes(x = Year, y = `Equivalent Area`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Equivalent Area (EA) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA))
ea_p
Compares the observed spatial density distribution to the density distribution expected from a homogeneously distributed population. Ranges from 0 to 1, with 0 indicating that the population was observed in only one ICES rectangle, and 1 indicating that the population density was uniformly distributed across ICES rectangles.
SPI <- spi(hlhh, yrs, qrs, species_aphia, stk_divs)[c(1,2,4)] # function needs tidying
head(SPI)
# A tibble: 6 × 3
# Groups: Year [6]
Year Quarter SPI.dur
<int> <chr> <dbl>
1 2000 1 0.435
2 2001 1 0.450
3 2002 1 0.455
4 2003 1 0.426
5 2004 1 0.500
6 2005 1 0.483
spi_p <- ggplot(data = SPI, aes(x = Year, y = `SPI.dur`)) +
geom_line() +
geom_point() +
theme_minimal() +
labs(title = "Spread of Participation Index (SPI) Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
theme(panel.border = element_rect(colour = "black", fill = NA)) +
ylab("Spread of Participation Index (SPI)")
spi_p
Plot all the the timeseries of indicators to visualise and compare trends.
# List of dfs
df_list <- list(cg, inert, eoo, ela, popr, poph, gni, D, sa, ea, SPI)
# Categorise indicators
loc <- c("CoG (x)", "CoG (y)")
ran <- c("Inertia", "EOO", "ELA")
occ <- c("POPR", "POPH")
agg <- c("Gini Index", "D95", "SA", "EA", "SPI")
# Merge all outputs into a single df
sidf <- Reduce(function(x, y) merge(x, y, all=TRUE), df_list) %>%
select(-c(nrects, nrects_p, no_haul.ids, pr_hauls)) %>% # remove some cols
rename(EOO = convex_hull_area,
POPR = PosAreaR,
POPH = PosAreaH,
ELA = `Ellipse Area`,
SPI = SPI.dur,
SA = `Spreading Area`,
EA = `Equivalent Area`) %>%
tidyr::pivot_longer(cols = 3:14, names_to = "Indicator", values_to = "Value") %>%
mutate(Type = case_when(
Indicator %in% loc ~ "Location",
Indicator %in% ran ~ "Range",
Indicator %in% occ ~ "Occupancy",
Indicator %in% agg ~ "Aggregation")) %>%
mutate(Indicator = factor(Indicator, levels = c(loc, ran, occ, agg)))
head(sidf)
# A tibble: 6 × 5
Year Quarter Indicator Value Type
<dbl> <chr> <fct> <dbl> <chr>
1 2000 1 CoG (x) 6.35 Location
2 2000 1 CoG (y) 55.8 Location
3 2000 1 Inertia 14.8 Range
4 2000 1 EOO 89.3 Range
5 2000 1 ELA 86.2 Range
6 2000 1 POPR 0.779 Occupancy
ggplot(data = sidf, aes(x = Year, y = Value)) +
geom_line(aes(colour = Type)) +
scale_x_continuous(breaks = seq(from = min(unique(sidf$Year)), to = max(unique(sidf$Year)), by = 2)) +
facet_wrap(vars(Indicator), scales = "free") +
theme_minimal() +
theme(
panel.border = element_rect(colour = "black", fill = NA),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6)
) +
labs(title = "Spatial Indicator Timeseries",
subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
ylab("Indicator Value")
R version 4.2.2 (2022-10-31 ucrt)
FLCore: 2.6.19
ggplotFL: 2.7.0
ggplot2: 3.4.4
Compiled: Tue Jan 23 15:46:15 2024